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FEASIBILITY STUDY OF THE NUMERICAL INTEGRATION 


OF SHELL EQUATIONS USING THE FIELD METHOD 
By Gerald A. Cohen 

Structures Research Associates, Laguna Beach, California 


SUMMARY 


The "field method" for the numerical solution of even order linear 
boundary value problems in ordinary differential equations is formulated. 

In essence, the field method converts the boundary value problem: into 
two successive initial value problems, each defined over the same one- 
dimensional domain as the boundary value problem. In the first of these 
problems the dependent variables are the "field functions", which appear 
in linear algebraic relations (the "field relations") satisfied by the 
dependent variables of the original boundary value problem. In general, 
for open branch domains the number of (scalar) field relations is equal 
to one-half the order of the boundary value problem. The field functions 
are obtained over the domain of the independent variable by a standard 
forward integration technique. In the second initial value problem the 
dependent variables are one-half of the original dependent variables, the 
remaining variables being given in terms of these by the field relations. 
The field functions appear as known terms in the differential equations 
of this problem, which is more properly termed a final value problem 
since it is solved by a backward integration over the domain of the 
independent variable. 

In this report the field method is developed for arbitrary open 
branch domains subjected to general linear boundary conditions. Although 
closed branches are within the scope of the method, they are not treated 
here. The numerical feasibility of the method has been demonstrated by 
implementing it in a computer program for the linear static analysis of 
open branch shells of revolution under asymmetric loads. For such problems 
the field method eliminates the well-known numerical problem of "long 
subintervals" associated with the rapid growth of extraneous, solutions. 
Also, the method appears to execute significantly faster than other 
numerical integration methods. 



INTRODUCTION 


Most problems in the response of axisymmetric shell structures may 
be reduced to the solution of equivalent linear statics problems. The 
earliest known application of numerical integration techniques to the 
solution of such problems was presented by Goldberg, et al. (ref. 1). 

In its most rudimentary form the mathematical problem is formulated as a 
two-point linear boundary value problem in ordinary differential equations. 
A set of linearly independent complementary solutions and a particular 
solution of the differential equations are obtained by a forward integra- 
tion scheme (such as Runge-Kutta) . Linear algebraic equations for super- 
position constants are then solved so that superposition of these auxiliary 
solutions gives the solution satisfying the boundary conditions. This 
method is, of course, nothing more than an attempt to use a well-known 
analytical approach for such problems as the basis of a numerical analysis. 

Although this approach works well when the auxiliary solutions are 
obtained in closed form, it is numerically ill-conditioned in the sense 
that for sufficiently long intervals of integration the superposition of 
the auxiliary solutions involves the small difference of large numbers 
at points remote from the initial point.* This problem of "long sub inter- 
vals" was subsequently circumvented by Cohen (ref. 2) and Kalnins (ref. 3), 
who divided the range of integration into suitably small subintervals. 

By using fresh initial values for the auxiliary solutions at the initial 
point of each sub interval, the growth of these functions is held to 
reasonable limits. Also, by allowing rings and other discontinuities to 
exist at the end points of each sub interval (i.e., the "boundaries"), the 
problem was simultaneously generalized to a multi-point boundary value 
problem. One difficulty with this approach is that the limit on sub- 
interval length depends on factors in addition to the properties of the 
structure itself, e.g. the circumferential wave number and, for eigenvalue 
problems, the eigenvalue shift. Also some problems may not be treatable 
simply because more subintervals are required than allowed by computer 
storage limitations. 

Later, Zarghamee and Robinson (ref. 4) further improved the method 
by proposing the use of initial values for the auxiliary solutions which 
imply satisfaction of the boundary conditions. Although this idea reduces 
the number of complementary solutions required by a factor of two, it 
does not eliminate the long sub interval problem. Their method was 
subsequently generalized to an arbitrary open branch domain [known as a 
"tree" in the terminology of geometric graph theory (ref. 5)] and general 
linear boundary conditions by Anderson, et al. (ref. 6). Cohen (ref. 7) 
further refined this method by improving the treatment of general boundary 


*If, instead, a numerical reintegration for the desired solution is 
attempted, the rapid growth of extraneous solutions associated with 
round-off errors will have the same effect. 
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conditions and generalizing it to domains which include a single closed 
branch. 

About the same time as the Zarghamee method was being developed , 
Jordan and Shelly (ref. 8) demonstrated a numerical integration technique 
for two-point boundary value problems which eliminates the long sub- 
interval problem. This method, which they termed the "field method", 
does not use complementary and particular solutions at all, but replaces 
the boundary value problem by two initial value problems to be solved in 
succession. The field method was later formulated by Miller (ref. 9) for 
two-point boundary value problems governed by a general ordinary linear 
differential equation of even order. 

The purpose of the present study is twofold: 1) to formulate the 

field method for a general even order multi-point boundary value problem 
defined on an arbitrary tree, and 2) to demonstrate the numerical 
feasibility of the method by implementing it in a computer program for 
static response of open branch shells of revolution subjected to arbitrary 
boundary conditions. 


SYMBOLS 


a,b,c,d 


p * p matrix coefficients of linear differential 
equations [eqs. (1)] 


B,D 


p x p matrix coefficients of linear boundary 
conditions [eq..-'. (2)] 


r (0) p (2) 

Cl 


meridional stretching and bending stiffnesses 


e 


p x p diagonal scaling matrix [eq. (14b)] 


f .»g 


inhomogeneous p x i matrices of linear differential 
equations [eqs. (1)] 


I 


p x p identity matrix 


L 

SL 


inhomogeneous p x 1 matrix of linear boundary 
conditions [eq. (2)] 

cylindrical length 


Mi 


meridional stress couple 


n 


circumferential harmonic number 


P,Q,S 


shell forces per unit of circumferential length in 
axial, radial, and circumferential directions 
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p 

R 

R 2 

r 

SC! ,SC 2 
s 


S , <j) , z 


t 


u,w 


X,y 


y.z 


A 

K 

C»h> v 


X 


Subscripts : 
0 


half-order of boundary value problem 
spherical radius of curvature 
circumferential radius of curvature 
small circle radius 

scale factors used as elements of e matrix 
arc distance 

meridional, circumferential, and normal coordinates 
of shell reference surface 

effective wall thickness 

p x p and p x 1 field function matrices [eq. (6)] 
axial and radial coordinates 

generalized force and displacement p x 1 response 
matrices 

net change across a vertex [eqs. (3) and (7)] 

B -1 D 

shell displacements in axial, radial, and circum- 
ferential directions 

meridional rotation 


value at initial shell edge 


Superscripts : 
T 


( )' 
( ) + 
( )' 
(~) 


matrix transpose 
d( )/ds 

value at vertex on exiting arc 


value at vertex on entering arc 
modified variable for singular arcs 



ANALYTICAL FORMULATION 


In order to formulate the field method in a general context, it will 
be convenient to introduce some elementary concepts from the theory of 
geometric graphs. Figure 1 shows the reference meridian of a hypothetical 
shell of revolution. The heavy dotted points depict boundaries of the 
shell, which are defined as one of the following types of points: 

(1) branch points 

(2) branch extremities 

(3) ring or ring load points 

(4) shell property or load discontinuity points 

The boundaries thus divide the meridian into a number of subintervals . 

In contrast to other numerical integration methods, no additional 
(artificial) boundaries are required simply to reduce subinterval length. 

A geometric graph (ref. 5) is defined as a set of points, called 
vertices, and a set of non-self-intersecting curves, called arcs, 
satisfying the following requirements : 

(1) Each closed arc contains precisely one vertex. 

(2) Each open arc contains precisely two vertices, viz. its end 
points. 

(3) The arcs have no common points, except for the vertices. 

It< is clear from figure 1 that if we identify the boundaries as vertices 
and the sub intervals as arcs, the reference meridian of a shell of 
revolution is nothing more than a geometric graph. 

The following graph terminology will be used: 

(1) Chain - a continuous sequence of arcs from an initial vertex 
to a terminal vertex. (A non-self -intersecting chain is said 
to be simple.) 

(2) Circuit - a chain whose initial and terminal vertices coincide. 

(3) Connected graph - a graph for which every pair of vertices is 
joined by at least one chain. 

(4) Tree - a connected graph which contains no circuits. 

In the present formulation of the field method, we shall confine 
our attention to boundary value problems defined over one-dimensional 
domains representable as (i.e., isomorphic to) a tree (fig. 2). It is 
important to recognize that every pair of distinct vertices of a tree 
is joined by precisely one chain, since connectivity implies the 
existence of at least one chain, whereas the absence of circuits implies 
the existence of at most one chain. It therefore follows that the cutting 
of any arc of a tree will disconnect the tree into two separate parts, a 
fact which is used in setting up the one -dimensional coordinate of the tree. 
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For each arc, the arc distance s will be used as the independent variable. 
Let us assume that the arcs have been ordered and oriented (with respect 
to the direction of increasing s) , but defer for the time being the manner 
in which this was done. 


Definition of Boundary Value Problem 

A system of ordinary differential equations of order 2p may always 
be written as a system of 2p first-order equations. If we group one-half 
of the dependent variables in the p x 1 matrix y and the other half in 
the p x 1 matrix z, the system of first-order equations may be written, 
for linear systems, as two matrix differential equations, viz. 

y' + ay + bz = f (la) 

z’ +; cy + dz = g (lb) 

where prime denotes differentiation with respect to s; a, b, c, and d are 
p x p matrix functions of s; and f and g are p x l matrix functions of s. 

Equations (1) are defined at every interior point of each arc of the 
tree. They are supplemented by linear boundary conditions, defined at 
the vertices of the tree, of the form 

BAy + Dz = L (2) 

where B and D are p x p matrices, L is a p x 1 matrix, and 

A y = l y + - l y" (3) 

Here, y + and y represent the values of y at the vertex on exiting (s 
increasing away from the vertex) and entering (s increasing towards the' 
-vertex) ..arcs respectively. As implied by the form of eq. (2), it is 
assumed that z is continuous at vertices.* A vertex and its boundary 
condition are said to be singular if the matrix B is singular; otherwise 
they are called regular. t 


v The majority of boundary value problems in mechanics are self- 
adjoint. This property is synonymous with being derivable from a 


*It is assumed on physical grounds that one-half of the dependent variables 
are continuous at vertices, and these shall be grouped in the vector z. 

The vector y may be viewed as a generalized "force" vector corresponding 
to the generalized "displacement" vector z, and -Ay represents the net 
external force entering the vertex. 

tA singular boundary condition implies a relationship between the 
components of the displacement vector z, i.e., kinematic constraint. 
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variational principle. Physically, this corresponds to systems which 
do not involve energy losses. For the system defined by eqs. (1) and 
(2), the conditions of self-adjointness may be shown to be 



(4a) 

(4b) 

(4c) 

(4d) 


where 


k = B -1 D 


(5) 


and the superscript T denotes matrix transpose. Although condition (4d) 
strictly applies only to regular vertices, it may be used also for 
singular vertices if a singular vertex is viewed as a limiting case of 
a sequence of regular vertices. Thus, for self-adjoint problems a 
singular vertex is the limit of a, sequence of regular vertices, for 
each of which eq. (4d) holds true.* 


The Field Method 

Since a cut at an interior point of any arc of a tree disconnects 
the tree into two separate parts, it is clear that the value of z at 
the cut can serve as a boundary condition which, along with the differen- 
tial equations (1) and boundary conditions (2) over one of the parts, 
determines y and z over that part independently of the corresponding 
data over the remaining part.i/ In particular, the value of z determines 
y at the cut, which in view of the linearity of the problem is expressed by 

y = uz + w (6) 

where u is a p * p matrix and w is a p * 1 matrix. Equation (6) is called 
a field relation since it is satisfied by the "field" of all possible 
solutions y,z (depending on the unused data on the remaining part). 


*Physically speaking, kinematic constraint may be approximated as close as 
one pleases by a sufficiently stiff spring. 

tit is assumed that giving z at the cut determines a unique solution, 
i.e., no nontrivial solution of the homogeneous forms of eqs. (1) and 
(2) satisfying the cut condition z = 0 exists. This is certainly the 
case if y = 0 for all solutions of the homogeneous forms of eqs. (1) and 
(2) subject to an v arbitrary homogeneous cut condition. Physically, this 
corresponds to a* situation where the only possible response under no 
load is a rigid -body displacement. 
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Correspondingly, u and w are called field functions. It follows that if 
we order and orient the arcs of the tree so that for every such. cut all 
of one part is described by s before any of the other part, then the 
determination of u and w will be an initial value problem. Such an. arc 
ordering and orientation will be achieved if at each interior vertex 
(i.e., one incident with two or more arcs): 

(1) there is precisely one exiting arc, 

(2) all entering arcs are ordered (i.e., described by ,s) before 
the exiting arc, and 

(3) the exiting arc is ordered immediately after the last entering 
arc. 

Implicit in these conditions is the requirement that one of the vertices 
of the first (last) arc is incident with no other arc, i.e., represents 
a branch extremity, and this arc is oriented away from (towards) this 
vertex. Assuming that the s-coordinate is set up in this manner, it is 
clear that, for a tree, eq. (3) reduces to 

Ay = y + - l y~ (7) 

since there will be at most one exiting arc for each vertex. 

Although u and w exist at interior points of each arc, this is not 
true at the initial vertex of an arc if the vertex is singular. (Such 
an arc itself will be called singular; otherwise an arc is regular.) 

The condition |b| = 0 is equivalent to the specification of a linear 
combination of the components of z at the vertex. This relationship 
plus eq. (6) would constitute initially on a singular arc p + 1 equations 
in the p unknown components of z. Compatibility of these equations 
implies a linear relationship between the components of y + . Since, 
however, y + depends on as yet unused data on the remaining part of the 
tree, it follows that eq. (1), i.e., u and w, cannot exist initially on 
singular arcs. From the foregoing, it is clear that singular arcs 
require special treatment. The discussion of singular arcs is postponed 
until after the next section, in which the basic method is presented. 

Trees without singular arcs .- Let us assume for the moment that all 
vertices are regular, except possibly one vertex corresponding to a 
branch extremity. In this case, the arcs can be ordered so that the 
sole arc incident with the singular vertex (if it exists) is the final 
arc, which is then oriented towards the vertex. Since the initial vertex 
of every arc is then regular, no singular arcs exist. This leads to the 
simplest form of the field method. 

In order to derive the differential equations for u and w, differen- 
tiate eq. (6) with respect to s, use eqs. (1) to eliminate y' and z', 
and eq. (6) to eliminate y to obtain 

(u' - ucu + au - ud + b)z + w' — ucw + aw + ug f~= 0 (8) 
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Since eq. (8) must be satisfied identically for all z (which depends on 
data at points of greater s, whereas u and w do not), it follows that 

u* - ucu + au - ud + b = 0 (9a) 

w' - ucw + aw + ug - f = 0 (9b) 

Equations (9) are the differential equations for u and w. 

The corresponding initial conditions for regular arcs are derived 


by substituting eq. (6) into eq. (2) to obtain 

(Au + B -1 D)z + Aw - B -1 L = 0 (10) 

where the synibol A has the same meaning as in eq. (7), viz. Au = u - | u 
and Aw = w + -! \ w“. Since eq. (10) is also an identity with respect to 
z, one obtains the initial values 

u + = -B -1 D + l u" (11a) 

w + = B -1 L + £ w" (lib) 

After integrating the initial value problem (9), (11) over the 
whole tree and storing the field functions u and w, eqs. (2) and (6) 
at the final vertex are solved simultaneously to give z there, viz. 

z = (D - Bu) _1 (L + Bw) (12) 


With this value as an initial condition, a backward integration of eqs. (lb) 
and (6), i.e. 

z' + (cu + d).z =i,g -* cw (13) 

is performed over the whole tree, using z-continuity at interior vertices. 
From the values of z so obtained, y is calculated from eq. (6), completing 
the solution. 

If the boundary value problem is self-adjoint, the numerical work 
is considerably reduced. In this case, the matrix u is symmetric so 
that only p(p + l)/2 of its component functions are independent.* This 
fact follows from eqs. (9a) and (11a) in view of eqs. (4). Thus in this 
common case, the total number of independent scalar differential equations 
contained in eqs. (9) and (13) is p(p + 5)/2, which compares with 
2p(p + 1) scalar first-order equations which must be integrated in the 
Zarghamee method. 


*In this case, u may be viewed as a '.'stiffness" matrix for that part of 
the tree, produced by a cut at s, which is fully described by smaller 
s-values. 
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Singular arcs .- Initially on singular arcs, the field relation (6) 
does not exist. Roughly speaking, one may say that the field functions 
u,w are infinite at such points. As discussed previously, the essential 
reason for this behavior is that a singular boundary condition implies 
a relationship between the components of z at the corresponding vertex. 

One is therefore motivated to make a transformation of variables to a 
modified z-vector. whose components are independent at the vertex. 
Consequently, on singular arcs new variables y and z are defined by* 

y = y (14a) 

z = z + ey (14b) 

where e is a constant diagonal p x p matrix required for dimensional 
homogeneity of eq. (14b). Substitution of eqs. (14) into eq. (6) shows 


that y,z satisfy the modified field relation 

y = uz + w (15) 

where u = (I + ue) -1 u (16a) 

w = (I + ue) _1 w (16b) 

Because of the symmetry of the transformation (14), eqs. (16) are inverted 
simply by replacing e by -e, i.e. 

u = (I - ue) -1 u (17a) 

w = (I - ue) _1 w (17b) 


In terms of these modified variables, eqs. (1) become 

y' + dy + bz = f (18a) 

2 ' + C? + d2 = g (18b) 

where 3 = a - be 
b - b 

C = c + e3 - de 

(19) 

d = d + eb 
f = f 

g = g + ef 


*An alternate transformation for singular arcs is given in Appendix A. 
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Equations (15) and (18) can be used to derive differential equations for 
u and w in exactly the same way as eqs. (1) and (6) were used to derive 
eqs. (9). Since eqs. (18) and (15) are identical in form to eqs. (1) 
and (6), respectively, the differential equations for u and w are in 
the same form as eqs". (9) with all variables replaced by corresponding 
tilde variables. 


The initial values of u and w on singular arcs may be obtained from 
eqs. (11) and (16). Considering the singular vertex as the limit of a 
sequence of nonsingular vertices, one has from eqs. (11) and (16) 


u + = lim [I + (£ u - B 1 D)e] 1 (£ u - B -1 D) 

I'B'bo 

= -B -1 D 

w + = lim [I + (£ u - B -1 D)e] -1 (£ w + B -1 L) 

I b bo 

= B -1 L 
where B = B - De 


(20a) 

(20b) 

(21a) 


D = D - B l u (21b) 

L - L + B £w~ (21c) 


The modified field functions u and w are calculated (and stored) 
on singular arcs by forward integration of eqs. (9) written for tilde 
variables, starting with the initial values (20). At the terminal vertex 
of singular arcs, it is convenient to replace u,w by u,w according to 
the reversion formulas (17). This is done so that the initial -values of 
the field functions for all arcs [eqs. (11) for regular arcs and eqs. (20) 
for singular arcs], as well as the terminal value of z [eq. (12)], may 
be computed the same way, regardless of whether preceding arcs are 
regular or singular. 

The backward integration on singular arcs is also done in terms of 
tilde variables, i.e., eq. (13) written for tilde variables is integrated. 
For this purpose, the value of z at the terminal vertex of a singular arc 
is replaced by z by first computing y there from eq. (6) (recall that u 
and w have been stored there) and then z from eq. (14b). After so doing, 
u and w at the terminal vertex are changed back to u and w [by using eqs. 
(16)], as required for the integration of the tilde form of eq. (13). 

From the values of . z obtained by the backward integration, y = y is 
computed from eq. (15) and z is computed from eq. (14b). 

For self-adjoint problems, it may be seen from eqs. (19) that the 
transformation (14) preserves the self-adjoint property of the differential 
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equations (1), i. e. , eqs. (4a-c) are satisfied in terms of tilde 
variables. Also, it is easily shown from eq. (16a). that the transforma- 
tion from u to u (and vice versa) preserves matrix symmetry. For eq. (16a) 
may be written, in the case of nonsingular u, as 

u = (u _1 + e) -1 (22) 

Since the inverse of a symmetric matrix is also symmetric, from eq. (22) 
symmetric u implies symmetric u. This result is .valid even if u is 
singular since symmetric singular u may be approximated as closely as 
one pleases by symmetric nonsingular matrices. In particular, since 
the value of u + [eq. (20a)] was derived as the limit of a sequence of 
symmetric matrices, u + is itself symmetric. As for u on regular arcs, 
integration of the differential equations for u preserves the symmetry 
of u . Hence, the simplifying conclusions drawn for self-adjoint problems 
on page 9 hold as well in the presence of singular arcs. 

A common type of boundary value problem originates as the minimization 
of a certain positive definite functional over a tree. An example of this 
type of problem is one-dimensional static response of an elastic structure. 
Such problems are self-adjoint, and the field relation (6) necessarily 
exists at all interior points of each arc, since the uniqueness condition 
mentioned in the footnote on page 7 is satisfied. To insure that the 
modified field relation (15) exists at all points of a singular arc 
(i.e., I + ue should be nonsingular), it is necessary to choose the 
diagonal matrix e [see eq. (14b)] positive definite, i.e., each of its 
nonzero elements should be positive. For such e, specification of z at 
a generic point s = s of a singular arc will uniquely determine y = y 
there, since this corresponds to cutting the tree at s and the attachment 
there of p stable elastic springs (to that part fully described by s < s). 
If e is not positive definite, at least one of the springs is unstable, 
which could lead to an instability (i.e., infinite u) for some s. 


SHELLS OF REVOLUTION 


A pilot computer program employing the field method to obtain the 
linear elastic response of open branch ring-stiffened shells of revolution 
subject to general harmonic mechanical and thermal loads has been written. 
For this class of problems, the tree over which the boundary value 
problem is defined represents the reference meridian of the shell. The 
differential equations (1) are eighth order so that the response matrices 
y and z are 4-element column vectors. The equations are ordered so that 
y and z are the force and displacement vectors (fig. 3) 

y T = r(P,Q,S,M 1 ) (23a) 

z T = (5>U»v,x) (23b) 
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where P,Q,S are forces per unit of circumferential length ref erred to 
fixed axial, radial, and circumferential coordinate directions x,y,4>, 
5,n,v are the corresponding displacement components. Mi is the meridional 
bending moment per unit of circumferential length, and x is the corre- 
sponding rotation. 

The matrices a,b,c,f and g [eqs. (1)3, as well as «..= B _1 D and 
B -1 L [eq. (2)] for ring boundaries, are given in Appendix B. Since this 
problem is self-adjoint, b,c and k are symmetric and d = -a^ [cf. eqs. 
(4)]. Corresponding shell and ring equations have been given previously 
(although not in precisely the same form) in reference 7. 

Solutions of several problems were obtained by the field method 
and compared to solutions obtained by the -Zarghamee method. The purpose 
of the numerical calculations was to: 

(1) uncover any practical problems in the implementation of the 
field method, 

(2) show that the "long subinterval" problem does not exist in 
the field method, and 

(3) compare the execution time of the field method with that of the 
Zarghamee method. 

Four basic shell configurations were studied: 

(1) clamped spherical cap (R/t = 91.4; n = 0) 

(2) branched conical shell (ro/t = 100; n = 1) 

(3) 140° sandwich cone (rg/t = 25.3; n = 0 and n = 1) 

(4) clamped-free cylinder (£/r = 1, r/t = 10; n = 0) 

Diagrams of the first three configurations are shown in figure 4. Here, 
n is the harmonic number considered, R is the spherical radius of curva- 
ture, r is a small circle radius, ro is the initial value of r, i is 
the cylindrical length, and t is the effective wall thickness (i.e., 
times the core depth for case 3, and the monocoque thickness for the 
other cases) . 

Each of these cases involve singular subintervals. In the first 
and fourth cases, this is due to the kinematic constraint at the initial 
edge. The second and third cases are free structures without rigid 
body constraint, but artificial rigid body constraint has been provided, 
in each case, at interior boundaries. The third case was previously 
used as the sample problem for the corresponding Zarghamee program 
(SRA 100) in reference 10, in which the variable pressure loading used 
is given. In the fourth case, the loading is uniform lateral pressure. 

During the course of the numerical evaluation of the field method, 
two practical problems were encountered. These are the choice of: 

1) the elements of the scaling matrix e, and 2) the storage locations 
for the field functions u and w. 
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Choice of the Scaling Matrix e 

In the general discussion of singular arcs, a scaling matrix e was 
introduced [eq. (14b)] without explicit definition, other than to say 
that it is a constant diagonal p * p (in this case, 4 x 4) matrix. In 
addition, it was deduced that for positive definite self-adjoint problems, 
to insure the existence of the modified field function u, all of the 
diagonal elements of e should be positive. Within this guideline there 
are an infinite variety of possible choices for e. Although the precise 
values of the diagonal elements of e have no effect on the solution, 
their order of magnitude determine the severity of the boundary layer in 
u following a singular boundary, and hence the rate at which the numerical 
integration will proceed. 

The first values tried for the diagonal elements e^ were based on 
the scale factors used in the supplemental initial conditions at singular 
boundaries in the Zarghamee method (ref. 7). These are eu = e ?2 = e 33 
= SCi - t/rC]/(P) and 644 = SC 2 - t/rCi ( 2 ) t where and Ci^ 2 ) are 

meridional stretching and bending stiffnesses, respectively. However, 
these values caused in sample problem 2 a severe boundary layer for u, 
and consequent slow forward integration, immediately following, the edge 
at point 7 (see fig. 4) when (artificial) rigid body constraints were 
imposed at that edge. After several empirical changes from these values, 
the values SCj - 10/Ci (0) and SC 2 - 3t/rCi( 2 ) we re chosen as optimal in 
the sense that for the sample problems studied they resulted in the 
least number of integration steps. The table below shows the number of 
derivative evaluations (four per Runge-Rutta step, whether accepted or 
not, plus one extra for each subinterval) during integration for three 
different sets of scale factors SCi and SC 2 . In each case, the upper 
value represents the forward integration and the lower value the backward 
integration. 


Scale factors 

Case 1 

Case 2 a 

Case 2^ 

Case 3 

Case 4 

SCi = t/rCi (O') 

— 

1441 





274 

SC 2 = t/rC^ 

— 

— 

— 

— 

234 

. lOOt/rCi (0) 

376 

425 

369 

421 

242 

t/rCx ( 2 ) 

312 

249 

217 

373 

194 

10/Ci ( o) 

344 

425 

361 

405 

234 

3t/rC! 

312 

249 

193 

373 

194 


Singular edge (point 7 of fig. 4) 

^Singular branch point (point 25 of fig. 4) 

The corresponding numbers for case 2 when the kinematic constraints are 
placed at the final edge (so that there are no singular sub intervals and 
the numerical process therefore does not involve the matrix e) are 377 
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and 177. Thus, it is seen that singular subintervals cause a small 
penalty when optimal scale factors are used. 

Not only is the numerical work reduced by the proper choice of 
scale factors for singular subintervals but at the same time an 
improvement in accuracy is generally realized. This is illustrated in 
the following table, in which the values of the dimensionless transverse 
shear stress resultant Q/pt and meridional stress couple Mi/pt 2 at the 
clamped edge of the clamped-free cylinder (case 4) are shown for the 
three field solutions of the preceding table, as well as for the 
Zarghamee solution. Here p represents the applied pressure load. 


Method 

Q/pt 


Mi/pt 2 


Field #1 

2.5181 

(3.16%) 

-2.8536 

(4.22%) 

Field #2 

2.4417 

(0.03%) 

-2.9738 

(0.18%) 

Field #3 

2.4413 

(0.02%) 

-2.9803 

(0.04%) 

Zarghamee 

2.4409 


-2.9792 



In this table the numbers in parentheses are the percent differences 
of the field solutions with the Zarghamee solution. Taking the Zarghamee 
solution as a standard,* the percent error in the field method is seen 
to be directly related to the number of integration steps required for 
the different scale factors. 

Although the third set of scale factors worked well in these test 
cases, no claim of universality is made for them. As programmed they 
are constant over the whole meridian, whereas it is only necessary that 
they be constant over each subinterval. One can anticipate that, for 
shells with large property variations along the meridian, it would be 
better to calculate different values of these-' scale factors for each 
singular subinterval. Also, for highly orthotropic shells for which the 
shear modulus is orders of magnitude different from the normal moduli, 
it may be desirable to introduce a third scale factor based on the shell 
wall shear stiffness. However, these refinements are not considered 
essential to demonstrating feasibility of the field method and are 
beyond the scope of the present study. 


Storage Locations for Field Functions 

When preparing the data deck for the computer program, the user 
must specify storage points for the response vectors y and z. The 
location of these points is arbitrary within the limitation that they 
are equally spaced within each sub interval. The backward integration 
for the z-vector is constrained to land on each of these specified points. 


*The Zarghamee solution is considered the most accurate of the four 
solutions since it does not require interpolation of calculated functions 
as does the field method, (in this regard, see next section.) 
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In the field method, it is. necessary to store the field functions 
u and w at the time they are calculated during the forward integration. 
Interpolation of these stored values is then made during the backward 
integration of eq. (13). In the initial version of the computer program, 
it was convenient to store these functions at the predesignated y,z 
storage points. Thus, in this version, the forward integration of eqs. 
(9) for u and w is also constrained to land on these points. However, 
because the variation of the field functions is of a totally different 
character than that of the response functions, this is generally a poor 
choice of storage locations for them. In general, the field functions 
have a narrow zone of rapid variation immediately following a boundary, 
but are otherwise slowly varying. In order to be able to interpolate 
accurately for intermediate values of u and w with the minimum number 
of storage points, the spacing’ of these points should vary with the 
rate of variation of u and w, i.e., the more rapid the variation of u 
and w, the closer together their storage points should be. Therefore 
the natural place to store them is at the end of each integration step, 
the size of which is automatically adjusted during execution according 
to their rate of variation. At the same time this would allow the 
forward integration to proceed at its own pace without being restricted 
by predesignated data points. 

In all of the cases studied, excellent agreement with the Zarghamee 
method was obtained for forces and displacements at the terminal shell 
edge. On the other hand, in some cases small errors of the order of one 
percent crept into the field method solution at other points due to the 
inadequacy of the specified y,z data points for the description of the 
field functions u- and w, resulting in interpolation errors for u and w 
during the backward integration. This is illustrated by the following 
table giving the dimensionless components of the z-vector obtained by 
each method at the initial (point 1 of fig. 4) and final (point 27 of 
fig. 4) edges of the branched shell (case 2, singular edge).* 


Edge 

Method 

5/t x 10 2 

n/t x io 2 

v/t x io 3 

X x io 3 

Initial 

Field 

-0.8352 

-1.4904 

0.3847 

-2.7215 


Zarghamee 

-0.8055 

-1.4595 

0.3833 

-2.6505 

Final 

Field 

-2.2100 

-2.8237 

4.2476 

3.1863 


Zarghamee 

-2.2096 

-2.8227 

4.2460 

3.1824 


This case showed the greatest loss of accuracy due to inadequate storage 
of the field functions of all cases studied. For example, the correspond- 
ing results for the n = 1 harmonic of the 140° sandwich cone (case 3), 
which has almost twice as many storage points (as well as fewer real 
boundaries) as case 2, are given in the following table. 


*This and all following comparisons are based on the same relative error 
control for the Runge-Kutta integration routine used in both methods. 
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Edge 

Method 

C/t * 10 2 . 

n/t x 10 2 

y/t * 10 2 

X x io 3 

Initial 

Field 

Zarghamee 

-0.9268 

-0.9353 

0.7079 

0.7110 

-0.0784 

-0.0786 

4.5021 

4.5386 

Final 

Field 

Zarghamee 

5.7869 

5.7874 

; -1.5461 
-1.5463 

1.4774 

1.4776 

0.5932 

0.5933 


Long Sub intervals and Execution Time 

In order to demonstrate that subinterval length has no effect on 
the numerical solution obtained by the field method, the clamped-free 
cylindrical shell configuration (case 4) was used. However, in order to 
put it into the range where several sub intervals would be required in 
the Zarghamee method (for which the axial'length of each sub interval 
should be less than approximately 5v / rt), the shell dimensions were 
changed so that £/r = 2 and r/t = 100. Also, for this case a comparison 
of the Zarghamee and field method integration times for the response 
due to uniforms fifth harmonic lateral pressure loading was made. For 
the Zarghamee setup, in order to avoid the long subinterval problem, the 
meridian was divided into four subintervals, each subinterval having 
four interior data points. For the field method only one interior 
boundary was used; this was specified at one-eighth the total length 
from the initial edge in order to change the spacing of data points 
qualitatively in accordance with the change in variation of the field 
functions. Four interior data points were specified in the first 
subinterval and six interior data points in the second subinterval. In 
order to avoid interaction with the choice of the scaling matrix e, 
singular arcs were avoided in this comparison by choosing the clamped 
edge to be the final edge, the initial edge then being free. Integration 
times (CPU times in seconds) for the CDC 6400 computer are shown in the 
table below. 


Integration 

Field 

Zarghamee 

Forward'd 

1.753 (290) 

5.156 

(516) 

Backward 

0.5411 (170) 

— 


Total 

2 .294 sec 

5.156 

sec 


The numbers given in parentheses .'are the number of derivative evaluations 
(four per Runge-Kutta step plus one extra for each subinterval) made 
during each integration. 

In spite of the rather meagre number of data points used to store 
the field functions, the accuracy of the field solution is good. The 
table below compares the dimensionless displacement components at the 
free edge obtained by each of the methods. 
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Method 

Kit. x .10 

n/t x 10" 1 2 3 4 

, v/t . 

X* 10 

Field 

4.4166 

2.2023 

-4.3587 

1.9081 

Zarghamee 

4.4222 

2.2033 

-4.3595 

1.9232 


CONCLUDING REMARKS 


The field method of solution of general even order linear boundary 
value problems defined on an arbitrary open branch one-dimensional 
domain has been formulated. The method has been implemented in a computer 
program for the static elastic response of open branch ring-stiffened 
shells of revolution subjected to general asymmetric (harmonic) loads. 

By studying specific sample problems the numerical feasibility of the 
method for axisymmetric shell problems has been demonstrated. For such 
problems the method eliminates the long subinterval problem of other 
numerical integration methods. In addition, it executes considerably 
faster than the Zarghameei method, which is the fastest numerical 
integration method currently in use. Numerical problems associated 
with kinematic constraints have been shown to be resolvable by the use 
of an "optimum" transformation of variables. It has also been shown 
that it is inadvisable to store the field functions at the same locations 
used to store the physical response functions. Rather, to obtain the 
greatest speed and accuracy for the least' storage, the field functions 
should be stored at the end of each integration step, the size of which 
is automatically adjusted during execution according to the rate at 
which they vary. 

Based on the demonstrated advantages of the field method over other 
numerical integration methods, it is recommended that the method be 
further developed with the object of obtaining a practical set of 
axisymmetric shell programs based on it. Such a set of programs should 
have the following features: 

(1) They should be applicable to shells with many closed branches. 

(2) They should be essentially free of numerical ill-conditioning 
problems (requiring, for example, double precision calculations). 
In particular, the user should not have to concern himself 

with questions of ill-conditioning when setting up a data deck. 

(3) Data deck preparation should be simple, at least no more 
difficult than that for other programs capable of treating the 
same problem. 

(4) They should execute faster and require less storage than 
existing programs for similar .‘problems. 
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APPENDIX A 


ALTERNATE VARIABLES FOR SINGULAR ARCS 

In place of eqs. (14) defining y,z on singular arcs, one could use 
the transformation 

y = y + ez (A-la) 

z = z (A- lb) 

where e is a constant diagonal p * p matrix required for dimensional 
homogeneity of eq. (A-la). In a similar fashion to the development of 
eqs. (15)- (21) for the tilde variables, the following analogous relations 
for the bar variables may be derived. The field relation is 

z = uy + w (A- 2) 

where u = (u + e) -1 (A-3a) 

w = -uw (A-3b) 

The inverse of eqs. (A-3) is 

u = u -1 - e (A-4a) 

w = -u -1 w (A-4b) 

The differential equations for y and z are 

y' + ay + bz = f (A-5a) 

z' + cy + dz = g (A-5b) 

where a = a + ec 

b = b - ae + ed 
c = c 

(A- 6) 

d = d - ce 

f = f + eg 

i = g 
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The differential equations for the modified field functions u and w are 


u 1 - ubu + du - ua + c = 0 (A-7a) 

w* - ubw + dw -K uf -r g = 0 (A-7b) 

The initial values of u and w on singular arcs are 

u+ = -D _1 B (A-8a) 

w + = D -1 L (A-8b) 

where B = B (A-9a) 

D = D - B(£ u _ + e) (A-9b) 

L = L + B l w~ (A-9c) 


With this modification for singular arcs, the calculation procedure 
is similar to that of the tilde modification (see pp. 10-12), except 
that in this case the backward integration on singular arcs is for y. 

The differential equation for y is obtained from eq.~; (A-5a) , which 
in view of the field relation (A-2) may be written as 

y' + (a + bu)y = f - bw (A-10) 

The integration of eq. (A-10) is started with the value of y, computed 
from eq. (A-la) , at the terminal vertex of the singular arc. 
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APPENDIX B 


MATRICES FOR AXISYMMETRIC SHELLS AND RINGS 


A 

E 


e x’ e y 


F x’ F y> F (j 


GJ 

I^j ly J Ixy 


k 


L 1 > l 2 


> ^f * 

N x ,N y ,N^ 


Xi,X 2 ,X 3 


0,(0), 0 2 (O) 

01 ( x ) ,0 2 (1) 


Symbols 

ring section area 

ring elastic modulus 

ring centroidal eccentricities 

harmonic amplitudes of ring force loads per unit 
of circumferential length 

ring torsional stiffness 

ring sectiori.-momehts- of ‘inertia 

ring stiffness matrix 

harmonic amplitudes of shell moment loads per unit 
of surface area in meridional and circumferential 
directions 

ring load vectors 

harmonic amplitudes of ring moment loads per unit 
of circumferential length 

harmonic amplitudes of shell force loads per unit 
of surface area in meridional, circumferential, 
and normal directions 

ring eccentricity matrix 

harmonic amplitudes of meridional and circumferential 
thermal force loads 

harmonic amplitudes of meridional and circumferential 
thermal moment loads 


0 


harmonic amplitude of ring free thermal strain 



orthotropic shell wall normal stiffness coefficients 
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orthotropic shell wall shear stiffness coefficients 
p ring centroidal radius 

Shell 

For elastic orthotropic shells of revolution subjected to symmetric 
n'th harmonic loads', the coefficient and load matrices of eqs. (1) are 
given below. For the definitions of the stiffness ( ^ ^ ) anc * l° ac * 
(X^,L^,0^( m )) variables in these matrices, see reference 7. 
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n 2 A 2 3 r, /R 2 n 2 X23r ' 2 /r n(r/R 2 -y 12 /r) n 2 A 2 4r,' /r 

-(Xi3+n 2 A 2 3/R2) ( r / R 2> “(*i 3 +n 2 A 23 /R 2 )r ' nr’ -(A 14 -fcn 2 A 2 4/R 2 ) 
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-rX 2 +(r/R 2 )L 1 +n[ (A 13 _+A 23 /R 2 )0 1 < 0 >-+-(A! 4+A 24 /R 2 ) 0j (l)-0 2 (O)-0 2 ( 1 ) /r 2 ] 
-rL 2 +r ' E A 2 3 © ! ^°^+A 24 0 1 



g » 


( r / R 2)[A 33 Q 1 (0)+x 3hQl (l-)j] 
r, fA 33 0j (O ) +A340i 

0 

A 34 0 iC o ) +Aii i +0j ( 1 ) 


Rings 

of"L bo «? arles > 


° f «n;r< r S: ;r fficleiit »- d ^ 


K e T ke 


7 ) 


where 


B' 1 L a -pT/„ 

(if + £ t - k £ e) 


e = 



o 2 (n 2 EI y « J)/p 2 n . 


e Wp z 


EA+n 4 EI / p 2 


Symmetric 


0 

1 


0 

0 


-e_ 


ne x/r ne y / r p/r 
0 


y' 

o 


n EI xy /p2 

n(EA-h 1 2 EI x /p2 ) 

n 2 (EA+ E r/ p 2 ) 


~n 2 (El y +Gj)/p] 

-" 2 EI xy /» 

-"*v 

EI y -Hi 2 Gj 


T 

% “ (pF y + nN nE » 

x “«y» pi* - n$ nE „ 

„ T y y x ’ pF <P ~ \> PNJ 

*t T - EA8(0,l,n, 0 ) 0 

£ e T = 6 (e x ,e y , 0 , 0 ) 
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